# 加载相关包
library(openxlsx)
library(reshape2)
library(tidyverse)
# 加载自编函数
source("./func/gee_run.r")
source("./func/summary_out.r")


# 读取数据，确定因变量、自变量、协变量的列
mice_data = read.xlsx("./output/0919_mice_data.xlsx", 4)
varsy = colnames(mice_data)[c(3:7)]        # 因变量
varsx = colnames(mice_data)[c(8:20)]       # 自变量
varsc = colnames(mice_data)[c(21:31)]      # 协变量

# 转换变量类型，注意：因变量为数值，自变量、协变量为因子
mice_data$id = as.character(mice_data$id)
mice_data$month = as.numeric(mice_data$month)
mice_data[varsy] = lapply(mice_data[varsy], as.numeric)
mice_data[varsx] = lapply(mice_data[varsx], as.factor)
mice_data[varsc] = lapply(mice_data[varsc], as.factor)
# 找出 8~24 月龄数据中填写次数至少两次的观测，步骤：melt_data => cast_data => 累加 8~24 的计数 => 提取计数大于 2 的 id 列表
extract_id = mice_data[, 1:2] %>%
             dcast(id ~ month) %>%
             mutate(sum_8_to_24 = `8` + `12` + `18` + `24`) %>%
             filter(sum_8_to_24 >= 2) %>%
             .$id
extract_data = mice_data %>% filter(id %in% extract_id & month > 6)
# 查看还剩下多少个个体
print(extract_data$id %>% unique() %>% length())
# View(extract_data)
# 将处理后的数据写入文件
write.xlsx(extract_data, "./output/0923_mice_data.xlsx", overwrite = T)


# 1. geeglm
## 1.1 删除自变量缺失的所有观测
extract_data = extract_data[complete.cases(extract_data[, varsx]), ]
 # 删除缺失值后，某些因子项的某个 level 计数变为 0 会导致 gee 报错；一种变通是重新计算因子的 level 来去除
extract_data[varsx] = lapply(extract_data[varsx], as.numeric)
extract_data[varsx] = lapply(extract_data[varsx], as.factor)


fit_table_1 = gee_run(
  extract_data,
  varsy = varsy,
  varsx = varsx,
  varsc = varsc[-5],
  file = "./output/0923_mice_10cv_del_birth_weight_c.xlsx"
)
fit_table_2 = gee_run(
  extract_data,
  varsy = varsy,
  varsx = varsx,
  varsc = varsc[-2],
  file = "./output/0923_mice_10cv_del_week_c.xlsx"
)

# give a fit_table, return a summary_plot_table
summary_out = function (
  fit_table
) {
  # summary_out
  ret = list()
  for (i in 1:length(fit_table[[1]])) {
    tb1 = list()
    for (j in 1:length(fit_table)) {
      tb1[[j]] = tbl_regression(fit_table[[j]][[i]], exponentiate = T, include = c(1))
      print(i)
    }
    ret[[i]] = tbl_merge(tb1, tab_spanner = varsy)
    
  }
  return (ret)
}

s1 = summary_out(fit_table_1)
s2 = summary_out(fit_table_2)
sect_properties = prop_section(
  type = "nextPage",
  page_size = page_size(),
  page_margins = page_mar()
)

save_as_docx(
  values = lapply(1, as_flex_table),
  path = "./output/0923_mice_10cv_del_week_c.docx", pr_section = sect_properties
)
save_as_docx(
  values = lapply(s2, as_flex_table),
  path = "./output/0923_mice_10cv_del_week_c.docx", pr_section = sect_properties
)
